Slow closure of denaturation bubbles in DNA: twist matters 
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The closure of long equilibrated denaturation bubbles in DNA is studied using Brownian dynamics 
simulations. A minimal mesoscopic model is used where the double- helix is made of two interacting 
bead-spring freely rotating strands, with a non-zero torsional modulus in the duplex state, = 200 
to 300 k^T. For DNAs of lengths = 40 to 100 base-pairs (bps) with a large initial bubble in 
their middle, long closure times of 0.1 to 100 /is are found. The bubble starts winding from both 
ends until it reaches a 10 bp met ast able state, due to the large elastic energy stored in the 
bubble. The final closure is limited by three competing mechanisms depending on and N\ arms 
diffusion until their alignment, bubble diffusion along the DNA until one end is reached, or local 
Kramers process (crossing over a torsional energy barrier). For clamped ends or long DNAs, the 
closure occurs via this latter temperature activated mechanism, yielding for the first time a good 
quantitative agreement with experiments. 



I. INTRODUCTION 

Since the discovery of the double helical DNA struc- 
ture by Watson and Crick in 1953 [1 , many studies have 
highlighted the role played by local DNA winding or un- 
winding in important biological processes, such as DNA 
replication, transcription or repair [2]. Biophysical exper- 
iments using single-molecule techniques |3] have shown 
that applying an external torque to DNA induces the 
formation of plectoneme or other structural changes 
Among them, the nucleation of a DNA denaturation 
bubble, a segment of several consecutive broken base 
pairs (bp), has been observed [5 and theoretically pre- 
dicted [BHIO] when a superhelical stress is imposed. 

In this paper, we focus on the role played by DNA 
torsional elasticity and twist dynamics in the sponta- 
neous closure of equilibrated large denaturation bubbles 
at room temperature. At first sight, once the bubble has 
been nucleated, for instance in vivo by the help of en- 
zymes, it should close almost instantaneously since the 
temperature is smaller than the denaturation one. How- 
ever, very large bubble lifetimes, in the 20 — 100 /is range 
for a 30 bps DNA, have been observed in in vitro exper- 
iments by Altan-Bonnet et at. [11 . These lifetimes are 
interpreted as closure times of the central bubble made 
of 18 AT (Adenosine and Thymine) nucleotides, flanked 
by two GC (Guanine and Cytosine) arms, known to be 
more stable. 

Several models p!2HT5] have studied the bubble breath- 
ing, i.e. intermittent and fast opening/closure of small 
bubbles, by considering an effective dynamics of the base- 
pairing states without focusing on the chain degrees of 
freedom. For instance, the Peyrard-Bishop model [16] 
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has been extended to consider twist degrees of free- 
dom pn4T9] . This model suffers two strong approxima- 
tions: (i) the helical axis is kept straight, i.e.. both bend- 
ing and the chain orientational entropy are neglected; and 
(ii) local breathing bubbles (or "breathers" ) emerge as lo- 
calized excitations of a non-linear wave equation which 
comes from a Hamiltonian dynamics with inertia [20] . 
However, since the dynamics of DNA in water is over- 
damped, these excitations have a lifetime of a few pi- 
coseconds [m |22] . These approaches are thus only valid 
at short time scale, such that the chain configuration can 
be considered as frozen, and cannot explain such large 
lifetimes as considered here. 

Other numerical works focus on the chain and inter- 
nal dynamics with various levels of coarse-graining, using 
molecular dynamics [23H26] or Langevin dynamics sim- 
ulations [27H3Q] . However they fail to capture the /is 
time scale time due to their high level of precision, or 
they are limited to short 30 bp long DNA [31]. In a re- 
cent paper [32j , we have proposed a simple coarse-grained 
model where two semi-flexible strands interact and form 
a planar "ladder" in the dsDNA form (without helicity). 
The coupling between base-pairing and bending elastic- 
ity was introduced through a varying persistence length 
equal to ^ds = 150 bp for dsDNA, and 4s = 3 bp for 
single-stranded (ss) DNA [33 . Closure times of 0.1 fis to 
4 fis^ following a scaling law of r ~ A^^-^, where N is the 
DNA length, were found, but still much smaller than the 
experimental closure lifetimes. 

In this paper, we improve this numerical model so that 
the two strands interwind to form a double helix in the 
double-stranded (ds) DNA state. The torsional modulus, 
Kcj)^ is chosen between 200 and 300 k^T in the dsDNA 
state (corresponding to torsional rigidity around 2.4 to 
4.5 X 10"^^ Jnm [3l|6l [34]) and taken to vanish in the 
bubble. We show that twist dynamics plays a key role in 
the closure of equilibrated large bubbles, which occurs in 
two steps. First, the large flexible bubble quickly winds 
from both ends {zipping regime), thus storing bending 
and torsional energy in the bubble, which stops when it 
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FIG. 1: (Color online) Snapshot of an equilibrated double- 
helix. The bending angle along each strand is ^o, po is 
the equilibrium base-pair distance, and h is the helical axis 
around which twist is defined. The imposed equilibrium twist 
between successive pairs is 00. 

reaches a size of ~ 10 bps. The closure of this metastable 
bubble depends on and A^: for low /t:^, an arms dif- 
fusion limited (ADL) regime is observed, as in previous 
ladder model, where the closure is controlled by the dif- 
fusive alignment of the two ds arms; for large /t:^ and 
not too large the bubble diffuses along the DNA and 
closes as soon as it reaches one DNA end [bubble diffusion 
limited (BDL) closure], with a closure time in r ~ N'^-^ 
for 40 < A/" < 100. For large /^^ and N or clamped ends, 
the closure is temperature activated (TA), which now ac- 
counts for the experimental observations |TT] . 

II. MODEL 

The DNA is modeled by two interacting bead-spring 
chains each made of TV beads (of radius 0.17 nm) of posi- 
tion r^. The Hamiltonian is H = Hg^^+H^f^+Htor+'^int, 
where the first two contributions are elastic energies of 
the strands j = 1,2 which include both stretching and 
bending energies 

N-l N-1 

The stretching modulus is Phz^ — 100, where = /cbTq 
is the thermal energy, Tq = 300 K is room tempera- 
ture, and ao = 0.357 nm. The bending modulus is 
large, ^hq = 600, to maintain the angle between two 
consecutive tangent vectors along each strand, 6>^, to the 
fixed value = 0.41 rad (see Fig[T]). Each strand is 
thus modeled as a Freely Rotating Chain (FRC) [36 . 
The third and fourth terms of H are the torsional energy 
and hydrogen-bonding interactions respectively. The tor- 
sional energy is modeled by an harmonic potential, 

N-l 

ntor=J2'^(^.-4>0f (2) 



where 0^ is defined as the angle between two consecu- 
tive base-pair vectors = r-"^^ — r-^^ and p^_^i (^o = 
0.62 rad). The stacking interaction between base-pairs is 
modeled through a hc^/y^i that depends on the distances be- 
tween complementary bases, hz^^^i = /^^[l — f {pi) f {pi j^i)] 
where f{pi) = [1 + erf (^^)]/2, and pi = |pj. Hence, 
= ^0 in the dsDNA state, and f<i(/y^i = in the ss- 
DNA one. We have chosen A' = 0.15 nm and = 1.5 nm 
and checked that a slight change in these values does not 
change significantly the results. The hydrogen-bonding 
interaction is modeled by a Morse potential, 

N 

^.„t = EA(e-^"^-2e-^) (3) 

where po = I nm, A = 0.2 nm, and f3A = 8 as in Ref. [32] . 

The evolution of Ti{t) is governed by the over-damped 
Langevin equation, integrated using an Euler's scheme, 

C^ = -^rM{rj})+m (4) 

where ( — Sirrja is the friction coefficient for each bead 
of diameter a with r] = 10~^ Pa.s the water viscosity. 
The random force of zero mean, ^^(t), mimics the ac- 
tion of the thermal heat bath and obeys the fluctuation- 
dissipation relation {^i{t) • Cj(^O) = Qk^TC,5ij 5{t — t'). 
Lengths and energies are made dimensionless in the units 
of a = 0.34 nm and k^T^ respectively. The dimension- 
less time step is 5t = 5tk^To/{a^Q^ set to 5 x 10~^ 
{St = 0.045 ps) for sufficient accuracy [32 . The equi- 
librium properties of this model DNA are described in 
the Appendix [Aj As an example, a typical equilibrium 
configuration of a 30 bp dsDNA is shown in Fig. [l] In 
particular, the fitted values for the dsDNA persistence 
length and the pitch are ^ds ^ 160 bp and p = 12 bp 
for Pf<i(iy = 300, which are comparable to the actual ds- 
DNA values (4s ^ 150 bp and p = 10.4 bp). The ss- 
DNA persistence length is 4s = 3.7 nm, compatible with 
experimental measurements [35j (see Appendix [A|. The 
initial bubble, of size L{t = 0) = N — 20, is created in 
the middle of the DNA by switching off Hint and then 
equilibrated for 3 ps. At t = the Morse potential is 
switched on in the bubble, and the dynamics is followed 
until the bubble closes. The cut-off value p* of the inter- 
base distance p defining closed (for p < p*) and open 
(for p > p*) states is fixed to 1.19 nm. Output values 
are then calculated every 1 ns, and samples are made of 
about 200 runs. Error bars are standard errors. 



III. BUBBLE CLOSURE DYNAMICS 

In Fig.[2|a) and (b) are shown typical evolutions of the 
bubble size, L{t)/L{0), for /3k.^ = 200 [Fig.|2];a)] and 300 
[Fig.[2|b)]. Two other geometric quantities related to the 
bending and twist stored in the bubble are shown: the 
scalar product rii • rie, where ni and rie are the tangent 
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FIG. 2: (Color online) Time evolution of the adimensional bubble size L(t)/L(0) (red), the scalar product between the two 
dsDNA arms axes, fii • fie, and the average twist angle per bp, Ac/), in the bubble for = 60 and (a) /SKcf) = 200; (b) /^/t:^ = 300. 
(c) Profile of the twist angles in the dsDNA just before (o) and just after the onset of the metastable regime (A), marked by 
arrows in (b) ((/)eq — 0.52 rad). 



vectors of both dsDNA arms (see snapshot in Fig.[4|, and 
the mean twist angle per base-pair inside the bubble 
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where io is the bubble first monomer. 

Two regimes can be clearly distinguished for any L{0) 
and /t:^: first a zipping regime, where L{t)/L{0) decreases 
rapidly, until it reaches a second metastable regime char- 
acterized by a stationnary L{t) = L ^ 10 bp. In the sim- 
ulations, we defined the onset of the metastable regime 
as the first time t such that L{t) = 11 bp. In the fast zip- 
ping regime, the initially flexible bubble closes due to the 
attraction between the two strands induced by the Morse 
potential. One example of the zipping dynamics is shown 
m Fig. [3]; a). The bubble size decays exponentially with a 
relaxation time on the order of 100 ns [102 and 208 ns for 
N = 70 and 100 respectively, see Fig. [sj^a)]. Indeed, dur- 
ing zipping, the two arms rotate in opposite directions 
in order to increase the twist of the whole chain and 
thus close base-pairs with an angular velocity, uj T/Co 
where T — 2A(j)Q^ 4/cBTrad is the driving torque and 
Co — 27r77Po^ — ^ kBTns is the rotational friction co- 
efficient of the arms {i is the arm initial length). We 
thus find u o:^ 1 rad/ns which induces zipping velocities 
V ^ puj/27r 2 bp/ns. This rough argument yields a 
consistent value with the zipping velocities measured in 
Fig. [sja) at short times. By defining the zipping time, 
Tzip, by L(rzip) = P = |[L(0)-Z], Fig.jsj^b) shows scaling 
laws, Tzip ^ with 1.4 < 7 < 1.5, as already observed 
with the ladder model [32 . Zipping occurs whatever the 
initial configuration, whether the two arms are aligned 
or not. 

The onset of the metastable bubble comes from the 
high 3D curvature of the two single strands inside the 
bubble, when its size reaches the ssDNA persistence 
length, Z ^ 4s- The two bubble single-strands are quite 
stiff at this scale. Either the arms are not aligned at the 
end of zipping and the elastic energy is both of bending 
and torsional nature, or they are aligned and it is only 
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FIG. 3: (Color online) (a) Semi-log plot of the bubble size 
L(t)/L(0) vs. time for 13k^ = 300, N = 70 (o), and N = 100 
(A). The solid lines are exponential fits, (b) Total zipping 
time as a function of P yielding an exponent of between 1.39 
(for I3k^ = 200), 1.49 (P^^ = 250), and 1.51 (P^^ = 300). 



of torsional nature. The non-zero twist at the onset of 
the metastable state (A0 ^ 0.2 to 0.3 rad) is created by 
the fast out-of-equilibrium dynamical closure of the bub- 
ble. This is illustrated in Fig. [2|c) showing the profile of 
the twist angle along the DNA just before (in red) and 
just after the onset of the metastable regime (in blue) for 
the simulation run shown in Fig. |2|b). It clearly shows 
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FIG. 4: (Color online) Evolution maps (open base-pairs are in black) for N = 60, and /^/t:^ = 200 (left), and /^/t:^ = 300 
(middle). Three snapshots of the DNA are shown, during zipping (bottom), in the metastable state (middle), and just before 
BDL closure (top). 



that the zipping stops as soon as the two domain walls, 
of approximately 5 bps width, "collide", increasing the 
twist angle value, and thus the twist energy, in the bub- 
ble centre. In brief, the zipping carries on as long as the 
elastic energy in the middle of the bubble is negligible. 
We have checked that the non-zero twist profile in the 
metastable state results from purely elastic properties of 
ssDNA (Appendix |B|. 

Depending on the value of Ai:^, the torsional contribu- 
tion of the elastic energy will be an energy barrier or 
not. Indeed, the closure mechanism, and therefore the 
dwell time in the metastable state, vary with For 

= 200 [Fig. [2|a)], A(f){t) increases smoothly until 
the bubble closes, whereas fii • rie increases from a nega- 
tive value, to a positive one in the metastable state. The 
bubble closure is thus controlled by the alignment of the 
two stiff arms, since closure occurs as soon as rii • rie — 1 
(ADL closure). This behavior has already been observed 
in the DNA ladder model [32 where no twist was present 
{hZ(f) = 0). The final closure was controlled by the rota- 
tional diffusion of one arm with respect to the other one: 
the metastable dwell time scaled with the DNA mean arm 
length, M, as r^^^ - with 2 < a < 2.4, and satu- 
rated at T]Pi^g for M > ^ds- We observe the same behav- 
ior for the helical DNA model with = 200, suggesting 
that, for this value, the twist does not play a significant 
role. As shown in Fig. [6]^b), we obtain r^^^^ ^ ^2.23 

= 200 (fitted solid line). The corresponding melt- 
ing map, shown in Fig.|4| illustrates that the bubble does 
not have sufficient time to diffuse far away from its initial 
position (the bubble diffusion coefficient along the DNA 
is i:> 1 bpVns). For I3i^^ = 300 [Fig. ^h)] however, 
the arms are almost aligned during the whole metastable 
state (but not necessarily during zipping). Moreover, the 
activation barrier to continue zipping being too high (see 
Eq. ([6| below), the fastest way to close is for the bubble 
to diffuse along the DNA until it reaches one end (Fig. [4| . 
This end opens to relax the torsion inside the bubble thus 
allowing a quick closure [43 . This BDL closure time is 
thus controlled by the one-dimensional diffusion along 
the DNA. We define precisely the arm alignment time by 
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FIG. 5: (Color online) Mean-squared displacement (in units 
of bp^) in the metastable regime of the bubble along the DNA 
vs. time for /^At^ = 300 and various N. 



the condition that rii -rie = 0.9. Figure |5] shows the mean- 
squared displacement of the bubble centre as a function 
of time for /^Ai:^ = 300. The bubble dynamics is purely 
diffusive, with a diffusion coefficient D 0.85 bp^/ns, 
almost independent of the DNA length. The final clo- 
sure is limited by the diffusion of the bubble towards 
one DNA end, leading to a dwell time in the metastable 
time, r^^t^ ^ {N/2f/{2D) ^ 0.15 A^^ ^^^^ ^^i^^ 
for I3i<i(f) = 200, in 25% of the simulation runs (50 over 
200) the bubble also closes using this mechanism (see 
Appendix [C]). The BDL regime starts to dominate for 
PhZ(f) > 230. Metastable dwell times, r^^^, and closure 
times, defined as the first time when the bubble closes 
completely, r^^^, are plotted in Fig. pi as a function of 
the DNA length N for Pi^^ = 300. The fit yields the scal- 
ing law Tmet ^ 0.06 A^^-^ which thus confirms the rough 
argument above (the prefactor changes due to a slightly 
different exponent). The total closure time follows the 
same scaling law r^^^ ~ N'^-^. We checked that, for 
l3K:(f) = 250, the exponent remains the same whereas the 
prefactor increases slightly to 0.075. 

A third type of closure exists: some trajectories show a 
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FIG. 6: BDL dwell (A) and closure times (o) {/3k^ = 300), 
and ADL dwell times (0) (/^/^^ = 200) with fits (see text). 
Horizontal arrows are the TA dwell times for 3 values of /Shicf). 
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Indeed, due to the connectivity of each DNA strand, all 
the base-pairs of such a small bubble can only close coop- 
eratively. To check this mechanism, we did simulations 
with clamped ends, which allowed us to avoid the BDL 
mechanism, and sufficiently large /t:^, to lower the ADL 
one. We clamped 10 bps on both DNA extremities (with 
Morse potential depth of 3^4/2) to represent either long 
or heteropolymer DNAs with GC rich sequences on each 
side, as in experiments [11]. For = 50 and fStz^j) = 300, 
out of 20 realizations, 12 of them did not close before 
100 /is and 8 of them closed in 52 /is on average. The 
bubble diffuses back and forth between the clamped arms 
several times and eventually closes. Figure clearly 
shows that the dwell time in this regime, i.e. the time 
actually spent by the bubble in the metastable state once 
the arms are aligned, follows an Arrhenius law 



^TA 



TO exp {EjkBT) 



(7) 



where tq is a prefactor almost independent of and the 
measured activation energy is ^Eg, O.IOa^^ — 17.6 [inset 
of Fig.[7|^b)]. By computing £^tor using Eq. (|6|, we find a 
comparable slope of 0.18. For /^/^^ < 200, the activation 
energy starts to saturate since we enter the ADL regime. 
Hence, from these simulations, it is clear that clamped 
DNAs, mimicking heteropolymer or long DNAs, take a 
long time to close, from tens to hundreds of jas. This is 
in quantitative agreement with the experimental results 
of Altan-Bonnet et al. [11 where an Arrhenius law was 
measured with Eg, = 7 kcal/mol ^ 11 /cbTq for = 30. 
Indeed extrapolating the inset of Fig. to this value 
yields a torsional modulus, f<i(/y = 280 /cb^o {C = 3.7 x 
10~^^ Jnm), a value consistent with observations [3j|34l 
l37] . Furthermore, the same activation energy value was 
measured in [11] for three different DNA constructs with 
an AT insert made of (i) a random sequence, (ii) a A track 
with its complementary T track, and (iii) a palindrome 
susceptible to form a cruciform. This is consistent with 
the scenario of a unique limiting step, that we show to 
be the formation of the 10 bps metastable bubble. 



IV. DISCUSSION 



FIG. 7: (Color online) (a) ADL and TA metastable dwell 
time distributions (free ends, /^k^ = 200, N = 100). Mean 
values are represented by vertical lines, (b) Arrhenius plot of 
the mean dwell time for clamped ends and pi^cf) = 220. Inset: 
activation energy vs. [N — 70). 



closure long after arms alignment but before the bubble 
reaches one end [Fig. [7|;a)]. This is a temperature ac- 
tivated (TA) closure, associated with the crossing of an 
activation energy barrier. Its torsional contribution, due 
to non-zero twist in the bubble, is: 



Eu 



K.4>,i) {(pi - 4>eq) 



(6) 



We performed several simulations for various N and 
/t:^, and constructed a "phase diagram", shown in 
Fig. (sja), representing the occurrence of the three clo- 
sure regimes in the (A^, a^^) plane. The methodology used 
to construct this diagram is described in Appendix [C| 
The definition of the various times contributing to the 
final closure time, Td, are sketched in Fig. [Sj with corre- 
sponding snapshots. For all the cases studied, the closure 
time, Tzip, is much smaller than the dwell time in the 
metastable state, r^et- In particular, as soon as the ini- 
tial bubble size, i^(0), is larger than L (and equilibrated), 
we expect the closure time to be essentially independent 
of L(0). 

The frontiers of the different regions (ADL, BDL, and 
TA) should be viewed as fuzzy since the diagram is estab- 
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FIG. 8: (Color online) (a) "Phase diagram" of the 3 mech- 
anisms of the bubble closure for a DNA with free ends (• 
corresponds to r^^^ or r^S^ = "^m^t; ♦ and ■ to simulations, 
see Appendix [Cj . For clamped (e.g. GC rich) ends, the BDL 
region is replaced by a TA one. (b) Sketch defining the three 
dwell times, r^^^, t^^^, and r^^t together with the zipping 
time, Tzip, and the closure time, Tci, and corresponding snap- 
shots. 



lished by comparing the metastable dwell times in each 
regime, r^^^^, r^^^, and r^^^, which are mean values 
of wide time distributions as shown in Fig. [7|^a). In the 
case of clamped (e.g. GC rich ^1 ) ends, the BDL re- 
gion merges into a TA one. For realistic DNA, one can 
assume /^/^^ > 200 [3^, 'F, '3T , which implies that only the 
BDL mechanism, for short DNAs, and the TA one, for 
long DNAs, might be observable. Furthermore, by ex- 
trapolating our results to very large both DNA with 
free ends and with clamped ends would have a bubble 
closure time which does not depend on N any more but 
is controlled by the local torsion, which provides a coher- 
ent picture of bubble closure for long DNAs inside the 
nucleus. 

A natural generalization of the model will be to con- 
sider the bubble sequence in the modeling, for instance 
by adjusting the parameter values in the interaction po- 
tential, Hint, and the torsional modulus profile, z^^^^, with 



the help of the Santa Lucia's nearest-neighbor model [38] . 
Taking into account the single-strand torsional elasticity 
would slightly increase the zipping time due to elastic 
resistance in the bubble but would not modify the occur- 
rence of the metastable regime. Finally, we did not con- 
sider hydrodynamic interactions in this work, and sup- 
pose the friction of the beads to be additive. The intro- 
duction of hydrodynamic interactions along one strand 
and between the two strands [39] might accelerate the 
closure, such as it decreases the relaxation times of sim- 
ple polymers. This work is in progress. 



Appendix A: Model DNA equilibrium properties 

This simple model captures most of the essential fea- 
tures of the system. The directionality is maintained by 
computing the sign of the determinant of (p^,p^^i,n) 
and then choosing the positive sign for a right handed 
helix (see Fig. [T]). The measured values of the geo- 
metric parameters are, after equilibration, = 1.20a, 
6>eq = 0.435 rad, and 0eq = 0.52 rad, that is slightly larger 
than the prescribed values ao = 1.05a, Oq = 0.41 rad, 
and (j)o = 0.62 rad due to thermal fluctuations and non- 
linear potentials entering the Hamiltonian. Moreover, 
our model DNA is a symmetric double-helix and not a 
double- helix with a major and a minor grooves. The 
ratio contour length/axis length is equal to 1.35 in our 
simulations, whereas it is equal to 1.7 for a real DNA [2]. 

The dsDNA persistence length, ^^s^ is computed using 
the method presented in Ref. for N = 150. The 
tangent-tangent correlation function C{s) = (t^+g • t^) is 
computed for each strand, where = t^/|ti| with = 
r^+i — is the unit vector connecting the two consecutive 
beads along a single strand. The correlation function is 
fitted, in Fig.|9|a), by the following theoretical expression 
(valid for a continuous helical chain) 



Cth{s) e 



(\ — u) cos 



/27r5 



(Al) 



where the persistence length -^p, the coefficient and 
the helical pitch p are fitting parameters. The fitted val- 
ues for the dsDNA persistence length and the pitch are 
^ds ^ 160 bp and p = 12 bp for p>K<^ = 300, which are 
comparable to the actual dsDNA values (£ds ^ 150 bp 
and p = 10.4 bp). Note that the equilibrium value of 
p is slightly larger than the prescribed one 27r/(/)o = 10. 
We have checked that the dsDNA persistence length is 
controlled both by bending and torsional potentials as 
they modify the local stiffness. For = 200, we find 
^ds ^ 100 bp. In the paper, we argue that the actual 
value for a real DNA is p>K^ = 280, yielding -^ds ^ 150 bp, 
as expected. 

We also estimated the persistence length of ssDNA, 4s ^ 
for N = 80. In fig.|9jb) is plotted the correlation function 
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FIG. 9: (Color online) (a) Correlation function C{s) for the 
simulated dsDNA {N = 150 bps). The solid line corresponds 
to Cth(s) given in Eq. (Al). (b) Semi-log plots of C{s) for 
the simulated ssDNA, allowing the determination of ^ss for 
three different values. (Inset) Fitted valu e fo r ^ss plotted 
together with the theoretical prediction Eq. ( A2 ) . 



C{s) for different values of in a log-linear plot. Due to 
the large value of the strength of the bending potential 
hZQ^ one can assume in a good approximation that 6> ^ 6>o, 
and £ss is purely controlled by the equilibrium bending 
angle Oq (freely rotating chain model). The correlation 
function is thus fitted by the exponential e~*/^^% which 
yields 4s = H bp, as shown in fig. |9jb). Moreover, we 
check that the ssDNA persistence length follows the law 



-a/ln(cosl9o) ^ 2a/0l 



(A2) 



for different values of [inset of fig. |9|b)], as ex- 
pected [36^. The value, 4s = 3.7 nm, is larger than 
the commonly accepted value of 1 nm. However the ss- 
DNA persistence length is not precisely known, since it 
has been shown experimentally [35 and theoretically [40 
that it varies with the salt concentration. Values of the 
order of 4 nm have even been found experimentally by gel 
electrophoresis [35]. The ssDNA persistence, 4s, cannot 
be modified in our numerical model without changing the 
pitch value because ^eq is a direct function of ^o- 



Appendix B: Geometry of the metastable bubble 

We have checked that the finite value of ~ 0.3, 
or the non-zero twist profile, in the metastable bubble 
results from purely elastic properties of ssDNA. We did 
some simulations to check the dependence of A^, in the 
metastable state, on the ssDNA elastic parameters and 
A^g. The procedure is as follows: we choose a snapshot of 
a metastable bubble for N = 60 and /^Ai:^ = 300. Then we 
switch the Morse potential off inside the 10 bp bubble and 
slightly decrease the temperature from T = Tq = 300 K 
to K (15 K every 30 ns). We observed that A0 remained 
equal to 0.3, thus confirming that the origin of this value 
is purely elastic in nature (and not entropic). 

Furthermore, to find the dependence on Acj) with Oq, 
we varied from 0.3 to 0.6 (without random forces for 
monomers belonging to the bubble). We found the lin- 
ear law A0 = 0.53^0 + 0.06. This is reminiscent of the 
3D bending of an elastic rod (see e.g. Ref. [41 ) with a 
spontaneous curvature Oq. 

Note that by decreasing the stretching modulus, a^^, 
we also observed an increase of A^ of 10% for /3hZs = 40, 
which might be a signature of the coupling between 
stretching and twisting as already mentioned in the lit- 
erature [42 . Finally, we have checked that by slightly 
changing the value of A' in the profile /^^^^ from 0.113 nm 
to 0.165 nm, we still observed the same metastable bub- 
ble size (data not shown). 

The 3D deformation of the single strands in the bub- 
ble comes form the constraint at their ends. They take a 
fluctuating helical conformation from which we can dis- 
tinguish two elastic contributions: (i) the bending is asso- 
ciated to the curvature of the central axis of the bubble, 
and (ii) the torsion is associated to the helical curvature 
of the strands, the central axis of the bubble remaining 
straight. Hence, the mean twist stored in the bubble in 
the metastable state results from a 3D bending of the 
bubble single-strands. 



Appendix C: Phase diagram construction 



In Figure p^a) are shown the dwell time distributions 
for BDL and TA closures. The procedure to measure 
them is as follows: for each trajectory, the ADL times and 
the TA or BDL times are measured. ADL times, r^^^, 
are elapsed times between the end of zipping {L{t) = 11) 
and the arms alignment (ni-rie = 0.9), BDL times, r^^^, 
are times between alignment and closure at one DNA 
end, and TA times, r^^, are times between alignment 
and closure inside the DNA (see Fig. [8|. One clearly 
observes the increase of the mean value and the spreading 
of the distribution with increasing N for the BDL case. 
In order to construct the phase diagram, we compare the 
average times of these distributions. All the data from 
simulations are given in Table [l| where hZ(f) is given in /cbTq 
and N in bps. For a given /^^ and the percentage of 
realizations belonging to BDL, TA and ADL are given. 
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N 


60 


70 


80 


90 


100 


BDL 


TA 


ADL 


BDL 


TA 


ADL 


BDL 


TA 


ADL 


BDL 


TA 


ADL 


BDL 


TA 


ADL 


200 


35.4 


49.4 


15.1 


28.4 


40.0 


31.5 


20.2 


49.2 


30.5 


12.5 


33.9 


50.5 


9.1 


31.7 


59.1 


210 


56.4 


35.9 


7.7 


44.4 


35.2 


20.4 


34.7 


43.5 


21.7 


28.3 


44.0 


27.7 


20.2 


45.0 


34.8 


220 


70.0 


22.7 


7.2 


56.8 


32.9 


10.1 


54.3 


37.5 


8.1 


42.6 


43.6 


13.7 


34.5 


38.0 


27.3 


240 


89.6 


9.2 


1.0 


89.6 


9.8 


0.5 


81.2 


16.2 


2.5 


74.8 


20.7 


4.3 


71.1 


23.9 


4.9 


250 


95.2 


4.7 


0.0 


90.5 


7.8 


1.5 


90.6 


6.7 


2.6 


90.1 


6.5 


3.2 


87.8 


9.4 


2.7 


300 


100 








100 








100 








100 








100 









TABLE I: Percentage of the bubble trajectories (DNA with free ends) following the closure mechanisms, ADL, BDL or TA. 



TA time ■ 
TV = 60 
N = 7Q I _-_ 
TV = 80 EI 
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FIG. 10: (Color onhne) (a) BDL dwell time, r^^^, distri- 
butions for various N and = 300, together with the TA 
time, T^^, distribution for free ends and /^a^^ = 200 [same 
as Fig. 4(a)]. (b) Evolution of the TA dwell time with 
for clamped ends. The values for /^a^^ = 200 and 210 (•) are 
extrapolated. 



The percentages are computed for ^ 200 realizations. To 
distinguish between ADL and TA closure mechanisms, 
we first calculated from the whole metastable trajectory, 

r^^L ^TA rj.^^^^ ^ADL ^ ^TA^ ^^^-^ ^^i^^ 

trajectory to belong to ADL case and vice versa. We used 
the same procedure to distinguish between BDL and TA 
closure mechanisms. 

The above data agrees with the phase diagram. We 
also did a few simulations for larger N\ 
• For = 200 & AT = 200 (♦ in phase diagram): 



56 realizations; 37 ADL closures, 1 BDL closure, and 
18 TA closures. This point is indeed slightly below the 
ADL/TA frontier line in the phase diagram of the article, 
as expected. The mean closure time is Td = 1.9 (±0.1) /is 
• For j3n^ = 250 & = 200 (■ in phase diagram): 56 
realizations; 1 ADL closure, 23 BDL closures, and 32 
TA closures. This point is thus almost at the frontier 
BDL/TA, as it can be checked in the phase diagram. 
The mean closure time is r^i = 4.7 (±0.4) /is. 

The few assumptions made in constructing the dia- 
gram are: ADL closure does not depend on ns] BDL 
times does not depend on [as shown in Fig. pi; and 
TA times are independent on N (local mechanism, see 
Fig. 4(b) of the main article). Furthermore, as one needs 
to know the TA times for pKff) = 200 and 210, since 
these values correspond to both ADL times and TA times 
(as observed in simulations), we did an extrapolation as 
shown in Fig[To]^b). All the times are plotted in the same 
Fig. [6] to extract the few data points for constructing the 
phase diagram. 

Relevant data points of the phase diagram are ex- 
tracted from Fig. [6j For example, the intersection be- 
tween a TA horizontal line and the ADL one for /^/^^ = 
200 gives = 52. Hence above = 52, the metastable 
bubble closes mainly through ADL. Likewise, the inter- 
section between a TA horizontal line for = 240 and 
the BDL one for pKff^ = 300 (we assume it is almost the 
same for 240) gives A^ = 120, which states that below 
N = 120 the bubble closes by BDL (most of the realiza- 
tions) and above which it closes by TA. 

Since TA times are given by = 
To ex.p{Ea{i^(f))/kBT) and assuming that ADL and 
BDL times do not depend on Tms = A^", equating 
both times yields the equation of the line separating 
TA and ADL or BDL regions in the phase diagram, 
/3hZ(f) = V -\- wlnN. By fitting the 5 data points for 
the frontier between BDL and TA regions, one obtains 
V = SS and w = 29. The fitted frontier line between 
ADL and TA (3 points) yields v = 101 and the same 
value for w. It is important to note that the frontier 
for low N between BDL, ADL, and TA is very fuzzy. 
Since for arm lengths larger than the dsDNA persistence 
length, M > £ds ^ 150, the ADL time does not depend 
on A^ any more [32 , the frontier becomes horizontal. 
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